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ABSTRACT 

Air flowing around the wing tips of an airplane forms horizontal tornado-like vortices that can be 
dangerous to following aircraft. The dynamics of such vortices, including ground and atmospheric 
effects, can be predicted by numerical simulation, allowing the safety and capacity of airports 
to be improved. In this paper, we introduce three-dimensional techniques for visualizing time- 
dependent, two-dimensional wake vortex computations, and the hazard strength of such vortices 
near the ground. We describe a vortex core tracing algorithm and a local tiling method to visualize 
the vortex evolution. The tiling method converts time-dependent, two-dimensional vortex cores into 
three-dimensional vortex tubes. Finally, a novel approach calculates the induced rolling moment 
on the following airplane at each grid point within a region near the vortex tubes and thus allows 
three-dimensional visualization of the hazard strength of the vortices. We also suggest ways of 
combining multiple visualization methods to present more information simultaneously. 
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1 Introduction 


Aircraft wakes represent potential hazards near airports. This hazard is so severe that it can control 
the required separation between aircraft, thus limit airport capacity. When multiple runways 
are used, crosswinds and density stratification can alter the trajectories and lifetimes of these 
wake vortices, producing hazardous flight conditions which can persist during subsequent flight 
operations. Thus improved forecasting techniques could result in safer airport operation and higher 
passenger throughput. 

In laboratories, experimental flow visualization techniques have been used to allow better un- 
derstanding of vortex interactions and of merging characteristics in multiple vortex wakes. For 
example, smoke and laser light sheet have been used to obtain both qualitative and quantitative 
information on the evolution of vortices. Up-to-date experimental methods remain the primary 
source of design information. On the other hand, because of the advances in computational tech- 
nology, numerical predictions of wake vortices have become feasible and are beginning to produce 
results consistent with experiments and physical observations. However, numerical predictions of 
wake vortex trajectories near the ground are still difficult. This is because the flow is unsteady and 
is characterized by at least one pair of moving viscous-trailing-line vortices, which respond to an 
essentially inviscid background, but are coupled to an unsteady, viscous ground-plane, boundary- 
layer region. The boundary-layer region can include separated flows during portions of the vortex 
wake history and thus create secondary vortices. 

During the development of more sophisticated numerical models, appropriate visualization 
methods are needed to monitor and verify the results from numerical simulations. Compared 
to the difficulties of experimental flow visualization, three-dimensional computer visualization is 
both flexible and repeatable. In this paper, we describe the visualization techniques that we have 
developed, along with development of numerical prediction models. 

Visualization of wake vortices can be divided into three steps: locating vortex cores, direct 
visualization of the vortices, and visualization of quantities associated with vortices, such as velocity, 
rolling moment, etc. The identification of a vortex core has been treated differently for different 
types of flows. A rigorous, widely-accepted definition of a vortex for unsteady, viscous flows is 
needed [9]. Direct visualization of vorticity fields is often misleading in recognizing the structure 
of vortex cores. 

Singer and Banks used both the vorticity and pressure field to trace the skeleton line of vortices 
in three-dimensional transitional flow [10]. Pressure field is used to correct the numerical errors 
that might be introduced during the integration of a vorticity line. A skeleton line passes the 
center of vortices. Pressure field is also used to determine the boundary of the vortex tube. Yates 
and Chapman examined both theoretically and computationally two definitions of vortex cores for 
steady flow: a minimum in the streamline curvature and a maximum in the normalized helicity 
[12]. However, the minimum streamline curvature does not give a viscous vortex core edge. That 
is, in our viscous flow calculations, the minimum streamline curvature only occurs at the center of 
the vortex, rather than at the edge of the vortex core. 

A somewhat vague definition of a vortex core is a vorticity-concentrated area characterizing 
the viscous aspect of the vortex. We are interested in the details of the vortex core such as its 
size and shape that change in time. Furthermore, for our time-dependent two-dimensional data, at 
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each time step, there is only one primary vortex and the center of vortex is well defined, where the 
maximum vorticity occurs. Our goal is to identify the boundary of the vortex core. 

We describe a tracing algorithm for identifying the vortex core at each time step using both 
the vorticity and velocity fields. Visualization of these vortex cores is done by tiling them into a 
tube-like surface, using a local surface construction algorithm. The surface is composed of Gouraud 
shaded triangle strips and can be efficiently displayed on a graphics workstation. It is colored by 
different scalar fields to examine different aspects of the vortex history. Two-dimensional techniques 
such as slicing and superimposed xy-plots are also used in conjunction with visualization of the 
vortex tube to present more information simultaneously. The techniques we have developed and 
the visualization results we have obtained help verify and understand the predicted flow field. 

An important quantity that may indicate vortex hazards is the induced rolling moment on 
the following airplane, which can be calculated using the unsteady two-dimensional predicted flow 
fields, combined with lifting line theory [6]. A novel approach presented here is to calculate the 
induced rolling moment on the following airplane at each grid point within a region near the vortex 
tube. The resulting data allow us to comprehend the hazard strength in three dimensions near the 
vortices with techniques such as direct volume rendering. Direct visualization of hazard zone in 
three dimensions can assist the design of aircrafts as well as the flight control at an airport. 

2 Numerical Model 

The difficulties related to numerical prediction of wake vortices near the ground have been discussed 
by Zheng and Ash [14]. We show briefly here the computational scheme for laminar cases. The 
coordinate system for this problem is shown in Figure 1. Because the problem is symmetric (without 
crosswind effects) in terms of the y-axis, only the first-quadrant of the flow field is computed. 
Treating the problem as an unsteady, two-dimensional flow, the velocity- vorticity formulation in 
Cartesian coordinates can be written as follows: 


where £ = 


dv du 

dx dy 


is 


, mn mn n , n, , 

dt dx dy dx 2 dy 2 

the vorticity. The stream function was governed by the Poisson equation: 


v 2 V> = -C, 

with 


( 2 ) 


IT ^ A 1 / 

U = — and V = . 

dy dx 


( 3 ) 


The vorticity transport equation is integrated by using an alternating direction implicit (ADI) 
scheme with upwind flux-splitting. The Poisson equation governing the stream function is computed 
using Swarztrauber and Sweet’s fast Poisson solver [11]. 

Computations were started when the two vortices were located at Xo = 1 , T 0 = 2 in non- 
dimensional units, where one unit is one half-span of the initial vortex pair. The initial vortex 
cores were assumed to have radii of 0.2, and symmetry was employed to reduce the size of the 
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Figure 1: Coordinate System. 


computational domain. A non-uniform grid, constructed by mapping a uniformly incremented 
exponential distribution into physical space, was employed, which was a 150x300 net. The grid 
points were packed in this manner to permit the vortex and boundary layer regions to be adequately 
resolved. The calculations were performed using the Cray II and the Cray YMP computer at NASA 
Langley Research Center. A typical laminar case, used in the following discussion, requires four 
megaword memory and four-hour CPU time on a Cray YMP to march 180 non-dimensional time 
units with At = 0.01. 

3 Vortex Cores 

The vortex core history is an important feature of wake vortex behavior. It is a way of checking the 
computational schemes used to calculate vortex problems and to characterize the vortex decay rate. 
In single vortex cases, as shown in Figure 2, the vortex core radius is defined as the distance from 
the vortex center to the point of maximum tangential velocity. In the current problem, a pair of 
vortices interact, and at the same time, interact with the ground boundary. The above definition of 
vortex core radius becomes vague and hard to apply in this case. The viscous interaction between 
the vortices and the ground causes deformation of vortex cores, and thus they are no longer circular. 
In addition, the complicated motion of the vortex in both vertical and lateral directions, and the 
use of Cartesian coordinates, makes the determination of the tangential velocity very difficult and 
somewhat misleading. 

We are investigating unsteady, two-dimensional laminar- vortex flows. The time dimension can 
be considered as the vortex axial dimension with constant speed motion in that direction. The 
unsteadiness, however, makes visualization of three-dimensional vortex cores difficult. First, the 
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vorticity values at the vortex core edge is different at each time step (see Figure 3), so that isovalue 
contours can not be used. Second, the same vorticity values may be generated near the ground 
boundary and thus even within a single time frame, isovalue contours may not reveal true vortex 
core shape. Finally, since the vortex core changes size and shape, the number of points included in 
the core changes. 

While it is natural to visualize the history of vortex cores by connecting them into a vortex 
tube, general-purpose visualization packages, like FAST [2], Tecplot [1], and others, cannot correctly 
construct the desirable vortex tubes. Some preprocessing is needed to extract the tube, after which 
we can perhaps make use of general-purpose visualization software. In the following sections, we 
describe a vortex core tracing algorithm and a surface construction method, which together form 
such a preprocessing step. 

3.1 Vortex Core Tracing 

The vortex core tracing algorithm developed here is based on the vorticity field, which is obtained 
directly with the vorticity-streamfunction formulation described in the previous section. The tan- 
gential velocity is now used only at the altitude of the vortex center (defined as the maximum 
vorticity point in the region of the vortex) to determine vorticity levels at the outer edge of vortex 
core. The tangential velocity at this altitude is simply the vertical component of the velocity field, 
which can be calculated using Equation 3, once the streamfunctions have been obtained. The left 
and right vortex core edges have thus been defined, where the maximum tangential velocity occurs. 
The average of the vorticity values at these two edge points is then used to establish the upper and 
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Figure 3: Vorticity Values at the boundary of each Vortex Core. 

lower boundaries of the vortex core at each lateral positions between the two edge points, marching 
from left to right. 

Figure 3 shows typical threshold values derived from one of our test data set for locating the 
boundary of the vortex core at each time frame. While the maximum vorticity is 7.93 and the 
minimum vorticity is -5.21 (near the ground), the threshold values vary between 0.052 and 1.792. 
Therefore, the use of a single isovalue would be inappropriate. This is also why we cannot use 
general-purpose visualization software. 

The tracing algorithm is called for each time frame (every 6 time units) saved for flow visualiza- 
tion, and therefore a time history of the vortex core can be shown. Figure 4 summarizes the tracing 
algorithm in C. Figure 5 shows a typical contour from data obtained at time = 180 when the sim- 
ulation ends. The points near the left and right edges are sparse due to the left-to-right traversal 
along the rectangular computational-grid that we used. The finer grid employed in y-direction is 
designed to resolve the ground boundary layer beneath the vortex. Redistributing points by using, 
for example, a cubic spline, would produce a better point set for the subsequent tiling task, probably 
resulting in a smoother surface. Nevertheless, this extra effort does not provide the computational 
scientist with additional information. Figure 5 also shows that at time = 180, the vortex core has 
grown in size and has been deformed, from a circle to an ellipse. 

3.2 From Vortex Cores to a Vortex Tube 

After the vortex core at each time frame is located, the next task is to tile consecutive vortex cores 
into a vortex tube. We can treat these consecutive vortex cores as a set of planar contours similar 
to those generated in medical imaging techniques such as computer tomography. The tiling problem 
for successive contours has been studied thoroughly and a summary of previous work can be found 
in the survey paper by Meyers et al. [8] on the subject of surfaces from contours. 

The efficiency of the surface construction process and the quality of surfaces generated is related 
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int center_x, center_y; /* center of a vortex core */ 

int left_edge_x, right. edge.x ; /* indices to edges of a core */ 
float vmin, vmax; /* most negative and positive velocity */ 

float threshold; /* for detecting core boundary */ 

for (t=0; tCnof.time.frames; t++) { 

F indCenterOf VortexCore (ftcenter.x , ftcenter.y) ; 

/* 

* find left and right edges of the vortex core 

*/ 

vmin=vmax=VerticalVelocity (center _x , center.y) ; 
for (int i=0; i<xdim; i++) { 

if (VerticalVelocity(i, center.y) < vmin) left.edge.x = i; 
else 

if (VerticalVelocity(i, center.y) > vmax) right edge x = i; 

> 

/* 

* find the boundary of the vortex core 
*/ 

threshold = (vorticity (left.edge.x, center.y) + 

vorticity(right_edge_x, center.y)) / 2; 

for (i=left_edge_x+l ; i<right_edge_x; i++) { 

/* find the upper boundary point */ 
for ( j=center_y ; j<ydim; j++) 

if (vorticity (i , j ) <=threshold) { 

AddUpperPoints(i , j ) ; 
break; 

> 

/* find the lower boundary point */ 
for (j=center_y; j>0; j — ) 

if (vorticity ( i , j ) <=threshold) { 

AddLowerPoints ( i , j ) ; 
break; 

> 

> 

> 


Figure 4: Vortex Core Tracing Algorithm. 
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Figure 5: Points Forming a Vortex Core. 

to the properties of the contours. Our contours (vortex cores) have the following features: 

1. Contours are parallel. 

2. Contours are usually well aligned. 

3. Contours are well shaped - no severe concavity. 

4. Contours have no holes. 

5. Each contour has a well defined center. 

6. The spacing between consecutive contours need not be exact. 

7. Each contour may have a different number of points. 

Feature 2 indicates that the contour (vortex core) changes its shape and position gradually in 
time. As to Feature 6, unlike, for example, in medical applications where the distance between 
consecutive contours needs to be accurate to make correct interpretation, the spacing can be arbi- 
trary as long as the user can determine the shape and size of the tube in the x and y dimensions. 
Depending on the speed of the airplane, the actual spacing can be very large with respect to the 
x and y dimensions. The total time units taken in the simulation may be equal to miles in the 
spatial dimension. So the spacing is selected for the ease of displaying and viewing with respect to 
the x and y dimensions. 
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Figure 6: Mesh of a Vortex Tube. 



Although Features 1-6 greatly simplify the tiling problem, Feature 7 introduces some difficulties. 
However, two existing metrics are reported in the literature for connecting consecutive contours. 
The first was introduced by Cook et al. [4], and is based on matching directions of points from each 
contour’s center of mass. The second metric was introduced by Christiansen and Sederberg [3] and 
is based on minimizing the length of diagonal between consecutive contours. Both of these metrics 
are Unear-time heuristic “greedy” methods that make the best choice locally, with no guarantee of 
finding a globally optimal result. In contrast, Fuchs et al. [5] proposes a globally optimal solution, 
which is computationally more expensive. In practice, Unear-time heuristic methods are usuaUy 
“good enough” and a lot faster when the contours have a lot of points. 

We have chosen the method proposed by Cook et al. and the result has been satisfactory. 
Figure 6 shows a typical mesh generated. The sparsity of points near the left and right edges of the 
vortex tube results in less desirable triangle combinations. However, this does not distract from 
the visualization of the history of the vortices. 

3*3 Visualization Results 

We appUed the above techniques to two-dimensional simulation results containing 31 time frames. 
The visuaUzation domain was a 140x280 non-uniform grid. The vortex tube generated was dis- 
played on a SiUcon Graphics Indigo 2. The three-dimensional graphics was implemented in GL 
and the user-interface was implemented in Motif. Four visuaUzation results are presented in Fig- 
ure 7. The upper-left image displays three color-contour sUces through the vortex tube shown in 
the upper-right image. The right-most sUce shows the vorticity field at the first time frame, which 
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touches the head of the tube. Color is mapped to vorticity and the colormap used is shown in 
the middle-left. The most negative vorticity value is mapped to blue, zero vorticity is mapped to 
white, and the most positive vorticity value is mapped to purple. Colors in between are linearly 
blended. The color of the vortex tube is mapped to vorticity values, £, in all the four images in 
Figure 7. Consequently, the head of the tube, in the right side of the image, is more red and purple 
indicating high vorticity values. The xy-plot superimposed on each color-contour slice shows the 
vertical velocity values along the lateral line passing through the center of the vortex core. The 
middle slice in the upper-left image shows that secondary negative vorticity has been induced by 
the primary vortex, after which separation occurs near the ground. 

The lower-left image shows the color-contour slice at the head of the vortex tube. The colors are 
mapped to kinetic energy and the colormap used is shown in the middle- right. Low kinetic energy 
values are mapped to white-blue and high to red-purple. Values in between are linearly mapped. 
Again, the xy-plot shows the vertical velocity. As you can see, the maximum kinetic energy region 
does not coincide with the vortex core (unlike in the single vortex case). The lower-right image 
displays the vortex tube from a different view, looking from the tail of the vortex tube. The tube 
becomes bigger as vorticity decreases with increasing time. 

In Figure 8, the upper image shows a volume rendered image of the vorticity field. Volume 
rendering offers direct visualization of complex structures and three-dimensional relationships. The 
image was generated by treating the stack of contours as a volume, and mapping vorticity values 
to color and opacity. By carefully selecting the color and opacity mapping, the volume rendered 
image also captures a vortex tube of similar shape. Unfortunately, this tube is less precise and is 
not identical to the one shown in Figure 7. The problem, which has been described in the previous 
section, is that a single mapping cannot capture the actual vortex tube. The true edge of the 
tube, where the maximum tangential velocity occurs, cannot be determined with direct volume 
visualization. The tangential velocity plot on the upper-right image in Figure 7 shows that edge. 

The lower two images in Figure 8 compare the vortex tube extracted previously by using the 
algorithm described in Section 3.1 with tubes extracted by using only an isovalue of vorticity. That 
is, each image shows as meshes the vortex tube superimposed with an isosurface of vorticity. The 
left image uses an isovalue of 0.6 and the right image an isovalue of 0.1. As expected, the isovalue 
tubes do not coincide with the vortex core tube. 

Further, compare the middle slice in the upper-left image in Figure 7 with the volume rendered 
image. The color contours display the vorticity field around the vortex core as well as near the 
ground plane. As mentioned before, the vorticity field generated on the ground may smear the 
displayed vortex core in the three-dimensional image. This phenomenon is shown in the volume 
rendered image too. The red vortex tube (positive vorticity values) is covered by a blue sheet 
(negative vorticity values) from the bottom. 

The vortex tube in Figure 7 gives an intuitive, three-dimensional display of the decay of vortex. 
This certainly resembles results obtained from experimental flow visualization and field observation 
(vapor trails). Additionally, computer visualization allows researchers to consider various atmo- 
sphere conditions such as turbulence, stratification and cross wind. The three-dimensional vortex 
tube will give a vivid picture of the behavior of the vortex under those different conditions. 
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Figure 7: Visualization of Vortex Tube. 
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Figure 8: Visualization of Vorticity Field and Vortex Tube. 
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4 Vortex Hazard - Rolling Moment 


In order to assess the vortex hazard, some measure of hazard strength is required. Since the 
computational domain is an unbounded quadrant, overall or global measures of circulation or 
velocity levels are of little value. The approach zone method used in [13] by Zheng and Ash could 
be meaningful in some sense, which calculated the circulation and kinetic energy in a bounded zone. 
However, the approach zone was defined somewhat arbitrarily, and the circulation and kinetic energy 
were not a direct measure of the vortex hazard. Therefore, a more systematic method has been 
developed by Zheng, Ash and Greene [15]. The rolling moment induced on following aircraft by 
the wake vortex flow field generated from leading aircraft is used as a measure of the hazard. 

In the rolling moment method, the down-wash velocity field, derived from Equation 3, is the 
computational result obtained from the numerical scheme in Section 2. Following aircraft are 
modeled as rectangular wings. In the current model, the influence of the following aircraft on the 
upstream flow field is neglected, as is the thickness of their wings. Then, employing PrandtPs lifting 
line theory [6], the induced rolling moment coefficient can be obtained using Fourier series expan- 
sions. Consider a spanwise circulation distribution on the following wing that can be represented 
by a Fourier sine series: 

N 

T(t) = 4-sWoo ^2 A n sin(n8 ), 

71 = 1 


where s is the half-span of the following airplane, is its forward speed and 6 is used to replace 
the physical spanwise coordinates of the following wing, x, defined as: 

cos9 = i. 

s $ 

As derived in [15], the rolling moment coefficient can be represented as: 

Ci = 7&A 2 
4 

where JSL represents the aspect ratio of the wing using the lifting line theory. We have 


1 r 

v(i} = 


1 (m dx . . ^ 


N 


1 x* — z 


sinQ 


YnAuSin{n9) 


n=l 


and if we can perform a Fourier sine expansion on u, which is obtained from the Navier-Stokes 
calculation using Equation 3, we can compute the rolling moment coefficient easily. The rolling 
moment data used in the following discussion required 1,200 CPU seconds on a Cray YMP. 


4.1 Visualization Results 

The example shown here is a case with the following airplane having the same wing span as that 
of the generating airplane. Since the flow field is symmetric (without crosswind effects), only half 
of the domain is displayed, as in the previous section. In the past, we calculated and displayed, 
with xy-plot, the history of rolling moment point by point. In order to visualize vortex hazard of 
the whole flow field, we use direct volume rendering of rolling moment coefficients, as depicted in 
the upper-left image of Figure 9. The colormap used is shown in the middle image in Figure 9. 
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The rolling moment coefficients visualized here, ranging from -0.189 to 0.054, are mapped to color 
and opacity. Negative values are mapped to green-blue, white represents no rolling moment, and 
positive values are mapped to yellow-red-purple. Higher absolute values are mapped to higher 
opacity. As a result, regions of lower rolling moment are more transparent and thus the region of 
hazard stands out. 

The upper-right image is a direct volume rendered image of positive-only rolling moment with 
the same view, from the head of the vortex tube. The lower-left image shows a different view, from 
the top (positive y). Note that, for all three images, only the region near the vortex tube is rendered. 
Use of the interactive volume visualization techniques described in [7], further exploration of the 
data can be made. It is seen that there are regions of positive rolling moment at both the left and 
right sides of the region of negative rolling moment. However, the left side is very small because it 
was annulled by the negative moment from the other side of the plane. 

It should be noted that the value at each point represents the rolling moment coefficient a 
following airplane experiences when the center of the wing span reaches that point in space. The 
following airplane is assumed to have the same size as the leading airplane. From our results, it is 
found that when the following airplane is near the region of the vortex core, the rolling moment is 
large (negative sign), and therefore the hazard is more significant. In fact, the location with largest 
negative rolling moment is the center of vortex. When it moves laterally away from the center, 
the rolling moment becomes smaller. If it moves laterally further to a position where only one side 
of the wing is influenced by the vortex, the rolling moment sign is changed, which is physically 
correct. Figure 10 shows such sign-change graphically. When the following airplane is in different 
positions relative to the center of the vortex, it experiences different direction of rolling motion. In 
Region 2, the following airplane moves into the region around the vortex core and counterclockwise 
rolling motion is induced on it. In Region 1, the port side of the following airplane experiences the 
rolling motion caused by the up-wash velocity of the vortex, and thus has a clockwise roll, while 
in Region 3 a clockwise rolling motion caused by the down-wash velocity exerts on the following 
airplane. 

It is observed that the maximum positive value is much smaller than the largest negative value, 
so one can conclude that the most hazardous region is the region around the vortex core. On the 
other hand, the kinetic energy contours in Figure 7 has lower values around the vortex center and 
higher values at both the left and right sides of the vortex core. Therefore, kinetic energy cannot 
be used to judge the hazardous regions. 

The lower-right image in Figure 9 displays two isosurfaces (blue and green) of rolling moment 
superimposed with the vortex tube (red) extracted previously. The isovalues chosen are 0.01 (green) 
and -0.01 (blue). As an example, the value of rolling moment coefficient -0.01 can be considered as a 
hypothetical hazard threshold for the following airplane. That is, the region inside the blue tube can 
be treated as where the induced rolling moment on the follower is beyond its controllable capacity 
and therefore the follower should avoid flying into this region. In addition, by superimposing tubes, 
the resulting visualization allows us to discern the relationship between the size of the vortex core 
and the hazard. It is easily seen that when the core size increases, the hazardous region decreases, 
which explains the importance of predicting correct core sizes in numerical simulations. 
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Figure 9: Visualization of Rolling Moment and Vortex Hazard. 
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5 Conclusions 

The study of the evolution of vortices behind a wing and its potential hazards is of considerable 
interest in aeronautic design and terminal control. In addition to experimental studies, numerical 
models have also been designed which might be used in practice for flight control at an airport. In 
this research, visualization techniques have been developed which allow the nature of numerically 
predicted wake vortices to be seen and analyzed. A vortex core tracing algorithm and a tiling 
method have been implemented which enable researchers to interactively examine the structure of 
vortex cores. Rolling moment, a quantity that may be used to measure the degree of hazard, is 
calculated throughout the domain of interest and visualized in three-dimensions using direct volume 
rendering. 

The comprehensive images generated by using these techniques suggest more intuitive ways of 
visualizing wake vortices and their corresponding hazard strength. These techniques, implemented 
with an interactive user interface, can be used not only for researchers to tune their numerical 
models, but also for flight control at airports. Consequently, our future work will include developing 
interactive visualization of vortex hazards and integrating the visualization techniques into the 
numerical simulation for real-time monitoring. 
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